{
//
// labFisInterfFit.C
//
// - open the Root file with data saved from interference experiment
// - get the microphone voltage vs. time experiment
// - create an empty monodimensional histogram
// - fill the new histogram with to voltages vs time from the first one
// - set the errors in the new histogram
// - fit a range of the histogram with a sine function
//

if (gSystem->AccessPathName("labFisInterfSave.root")) {
  fprintf(stderr, "Cannot find file labFisInterfSave.root!\n");
  return;
}
f = new TFile("labFisInterfSave.root");

hsMicroTimeProf = (TProfile *)f->Get("hsMicroTimeProf");

UInt_t bins = hsMicroTimeProf->GetNbinsX();
Double_t xlow = hsMicroTimeProf->GetXaxis()->GetXmin();
Double_t xupp = hsMicroTimeProf->GetXaxis()->GetXmax();

hs = new TH1F("hs", "Microphone voltage vs time", bins, xlow, xupp);

//
// errors
//
// here we are fitting a series of ADC measurements
//
// if there is no noise the error of one ADC measurement can be set
// as 1/sqrt(12) ADC counts (the sigma of a uniform distribution between
// x and x+1
//
// if there is noise, one can repeat the the ADC measurement of a constant
// voltage and measure the sigma of the ADC count distribution: this
// is a good estimate of the error to be assigned
//
// Here we assume that the error is 4 ADC counts
// For the interference experiment the ADC range is -5V .. + 5V
// and the ADC counts are from 0 to 0xffff, according to this the
// ADC count error is translated to voltage error
//

Double_t error = 4. * 10. / 0xffff;

//--- build histogram to fit, with data points and errors
for(UInt_t i=1; i< bins+1; i++) {
  //--- set bin i content (BEWARE: does not correspond to filling)
  hs->SetBinContent(i, hsMicroTimeProf->GetBinContent(i));
  hs->SetBinError(i, error);
}

//--- fit between 0.00425s and 0.006s
Double_t fitLow = 0.00425;
Double_t fitUpp = 0.00600;
Double_t fitUpp = 0.00435;

// 
// fit function = offset + A * sin( 2*pi*freq * x + phi )
// parameter offset accounts for constant voltage offset
//

TF1 *fun = new TF1 ("fun", "[0]+[1]*sin(atan(1)*8*[2]*x+[3])", fitLow, fitUpp);
fun->SetParNames( "offset","A", "freq", "#phi");
//--- initialize parameters to sensible values
fun->SetParameter( 0, 0);
fun->SetParameter( 1, 0.05);
fun->SetParameter( 2, 40000);
fun->SetParameter( 3, 0);

//
// set proper limits of parameter A,
// otherwise the fit does not converge well
//
// one could scan the data beforehand to get proper limits
//

//
// since the fit has trouble converging, a preliminary fit is done
// on just few points, in the range 0.0425s - 0.0435s
//

fun->SetRange(0.00425, 0.00435);

hs->Fit("fun", "r");

//
// the above preliminary fit converges well and finds good starting point
// for the parameters. The fit is repeated on the extended range.
//

fun->SetRange(0.00425, 0.00600);

hs->Fit("fun", "r");

gStyle->SetOptFit(111);
hs->Draw();

printf("\n");
printf("Fit results for f(x) = offset + A * sin( 2*pi*freq * x + phi\n");
printf("chi2            %6.4f for %d degrees of freedom\n",
       fun->GetChisquare(), fun->GetNDF());
printf("constant offset %10.4f +- %6.4f mV\n",
       1000*fun->GetParameter(0), 1000*fun->GetParError(0));
printf("amplitude       %10.4f +- %6.4f mV\n",
       1000*fun->GetParameter(1), 1000*fun->GetParError(1));
printf("frequency       %10.4f +- %6.4f Hz\n",
       fun->GetParameter(2), fun->GetParError(2));
printf("phase           %10.4f +- %6.4f\n",
       fun->GetParameter(3), fun->GetParError(3));

}
